'Trig  FILE  CGFV  AD  A 1  0  1 


\ 


AN  ADAPTIVE  ABMA  SPECTRAL  ESTIMATOR:  PART  2 

Jaoas  A.  Cadzow  and  Randolph  L.  Moses 
Department  of  Electrical  and  Computer  Engineering 
Arizona  State  University 
Tempe.  Arizona  85281 


iViV- 


Abstract  -  This  paper  presents  a  recursive  algorithmic 
Implementation  of  the  prewindowed  high  performance 
method  of  ABMA  spectral  modeling  as  described  In  Parc  I. 
This  algorithm  provides  updates  of  the  ARKA  models 
optimal  autoregressive  coefficients  (In  actuality  pre¬ 
diction  errors)  as  each  new  data  point  becomes 
available.  The  algorithm  Is  computationally  fast  in 
the  sense  that  It  requires  0(p)  multiplications  and 
additions  for  each  update.  It  Is  shown  chat  this  fast 
recursive  algorithm  may  be  Implemented  using  a  lattice 
filter  arrangement,  and  It  therefore  exhibits  several 
of  the  "nice"  properties  associated  with  lattice  type 
tlgorlchms  such  as  numerical  robustness  and  good  con- 
/ergence  properties. 

I.  INTRODUCTION 

In  Parc  1  of  this  paper  we  described  an  algorithm  for 
■btalnlng  an  estimate  of  Che  power  spectral  density 
ssoclated  with  a  given  time  series  {x(n)f.  In  partl- 
ular.  the  following  ABMA  (p,q)  spectral  density  model 
as  hypothesized 
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A  closed  form  algorithm  for  estimating  this  taodel's 
a^  and  bj  parameters  was  Chen  developed  in  which  the 
finite  sec  of  time  series  observations 


x(l),  x(2).. 


x(n) 


(2) 


were  used  in  Che  parameter  selection  process.  This  so 
called  "high  performance"  algorithm  operates  on  a  data 
block  of  length  n  to  obtain  the  model's  coefficients 
In  a  single  computational  effort.  It  Is  therefore 
called  a  block  processing  algorithm. 

There  are  many  situations,  however.  In  which  a  block 
processing  algorithm  for  spectral  estimation  Is  not  an 
appropriate  cool.  In  a  variety  of  applications,  Che  data 
measurement  Is  an  ongoing  process  and  it  is  therefore 
desirable  to  recursively  update  Che  autoregressive  and 
moving  average  parameters  as  each  new  data  point 
becomes  available.  This  capability  is  of  particular 
Importance  in  chose  cases  where  one  wishes  to  adaptive¬ 
ly  model  the  spectrum  of  a  long,  ongoing  time  series. 
Algorithms  with  this  recursive  updating  capability  are 
called  "recursive  algorithms". 

Recently,  several  fast  recursive  spectral  estimation 
algorithms  have  been  developed  [1]-[71.  Host  of  these 
algorithms  are  based  on  so  called  least-squares  AR 
spectral  estimation  methods  [8],  although  a  few  ABMA 
algorithms  have  also  been  developed  [4]  i  [6],  These 
slgorlchms  seek  to  minimize  a  prediction  error  vector 
In  order  to  obtain  the  desired  spectral  esclmstes. 

In  this  paper  we  develop  a  more  affective  ARMA 
algorithm  that  efficiently  computes  the  optimal  auto- 
ragrasslve  coefficients  by  recursively  updating  a  set 
of  prediction  error  elements  as  each  new  data  point  la 
obaerved.  This  recursive  algorithm  is  based  on  the 
prevlndowed  high  performance  method  as  described  In 
Pert  1,  and,  therefore  Is  predicated  on  the  approxi¬ 
mation  of  the  ARMA  model's  underlying  Yule-Walker 
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equations.  Although  similar  algorithms  may  be  derived 
for  the  unmodified  and  for  the  two  other  modified 
versions  of  the  high  performance  method,  the  recursive 
algorithm  based  on  the  prewindow  version  Is  a  bit 
easier  to  derive  and  Is  characterized  by  a  "fast  start¬ 
up"  capability  In  that  spectral  estimates  are  possible 
with  as  few  as  two  data  points.  The  recursive  algoritln 
herein  presented  is  computationally  efficient  In  the 
sense  that  0(p)  multiplications  and  additions  are 
required  to  update  the  "necessary"  parameters  as  each 
new  data  point  Is  observed.  This  paper's  recursive 
algorithm  was  originally  developed  in  [91.  A  more 
straightforward  derivation  is  herein  presented  which 
provides  a  greater  degree  of  Insight.  In  addition,  a 
lattice  filter  Implementation  of  this  algorithm  is 
developed.  Moreover,  because  of  this  ladder-type 
implementation,  this  algorithm  is  characterized  by 
several  other  nice  properties  associated  with  ladder 
algorithms  such  as  numerical  robustness  and  good 
convergence  properties  [10] ,  [11]. 

II.  THE  PREDICTION  ERROR  VECTORS 

The  recursive  update  equations  herein  presented  do 
not  explicitly  update  the  ARMA  model's  autoregressive 
coefficients  in  obtaining  optimal  updated  spectral 
estimates.  Instead,  a  set  of  "equivalent"  parameters 
known  as  prediction  errors  are  updated.  In  this 
section,  we  discuss  the  relationship  between  the  pre¬ 
diction  errors  and  autoregressive  coefficients. 

As  outlined  In  Part  1  of  this  paper,  the  optimal  p'*’ 
order  set  of  autoregressive  coefficients  for  the  pre- 
wlndomd  version  of  the  high  performance  method  are 
obtained  by  solving  the  following  system  of  p  linear 
equations  In  p  unknowns 


[Y^  X  la  +  Y'*’  x  •  8 

n*p  n»p  —p  n,p 

(3) 

where 

X  •  x(2),.  .  x(a>]' 

(4a) 

(4b) 

x„.p  -  ... 

(4c) 
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In  which  6  is  the  zero  vector  and  S  Is  the  do%m  shift 
operator. r  Here  the  dagger  symbol  (t)  denotes  complex 
conjugate  transposition  and  the  prime  symbol  (')  denotes 
transposition.  The  subscripts  p  and  n  explicitly  ,  — 
Indicate  that  the  denominator  order  of  the  spectral 
model  of  equation  (1)  is  p  and  chat  n  data  points  are 
available.  Whenever  this  explicit  information  la  not 
needed,  we  will  uee  x,  X,  and  Y  In  place  of  Xj,, 

Xn  p,  and  Vp  a,  respectively.  [ 

it  Is  recalled  from  Part  1,  chat  Che  high  performance 


ARMA  modeling  approach  is  predicated  on  approximating  t 

^In  this  paper,  matrices  are  denoted  by  capital  letters 
(e.g.  X),  vectors  are  denoted  by  underlined  lower  case 
English  letters  (e.g.  x)  and  scalars  are  denoted  by 
lower  case  Greek  letters  (e.g.  o) .  Moreover,  the  down 
shift  operator  S  is  defined  by 
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Yule-Walkar  aquation*  whare  c>p.  Upon  axaalnacion  of 
axprasaton  (3) ,  It  Is  apparanc  chat  wa  hava  hara  ra- 
scrlccad  c>p.  Thla  raacrlccion  Is  raqulrad  In  ordar 
to  facllltaca  ch*  davalopaanc  of  ch*  fast  rscurslva 
algorithm.  Unfotcunacaly,  by  caqulclng  tap,  cha 
■pactral  aaclMClon  parfonanca  suffars  in  comparison 
to  that  achlavad  with  largsr  valusa  of  c.  as  tha  data 
langch  n  Incraasas,  howavar,  this  parformanca  dagrad- 
aclon  dlmlnlahas  and  typically  Is  of  an  Insignificant 
nacura.  Thla  la  Indaad  fortunaca  slnca  It  la  praclsaly 
for  long  data  langth  casas  that  tha  racuralve  algorithm 
would  most  llkaly  be  ucUlzad. 

Under  the  assumption  chat  T+X  Is  nonsingular,  ch* 
optimal  autoregressive  coefficient  vector  which  satis¬ 
fies  expression  (3)  Is  given  by 

a  •  -  -(Yt  „]'V  X  (5) 

— p  '  n,p  n.p-'  n,p-n 

In  what  Is  to  follow.  It  Is  beneficial  to  Interpret 
this  autocorrelation  coefficient  salecclon  procedure 
from  a  prediction  error  viewpoint.  Namaly,  wa  may  ra- 
fomulace  expreealon  (3)  as 

Y+  f*  -  8  (6) 

n.p  -p,n  - 

in  which  „  Is  the  so-called  forward  prediction  error 
as  specified  by 


It  is  referred  to  as  Che  forward  prediction  error  since 
Its  component  can  be  interpreted  as  being  tha  error 
resulting  from  a  prediction  of  the  element  x(k)  by  a 
linear  combination  of  the  p  most  recent  time  series 
elements  x(k-l) ,  x(k-2) , .  .  . ,  x(k-p) . 

The  optimal  autoregressive  coefficient  vector  (5)  can 
be  then  associated  with  an  auxiliary  minimization 
problem  involving  the  prediction  error  vector.  Nawsly, 
It  is  readily  shown  chat  this  optimal  vector  minimizes 
the  following  quadratic  functional 

g(a  )  -  f?  „"W  f*  „  (8) 

— p  •  p » n  ““P  *  n 

where  W  Is  the  nxn  positive  semideflnlte  matrix 
specified  by 

W  -  Y  Y  ^  (9) 

n.p  n.p 

To  reinforce  this  prediction  error  Inccrprecaclon,  let 
us  define  cha  following  estimate  of  vector  ^ 

X  -  -X  a  (10) 

n  Rg p 

which  in  turn  generates  the  forward  predicting  error 
■*p « n 

Upon  substitution  of  expression  (3)  Into  (10)  the  opti¬ 
mum  forward  prediction  error  vector  la  given  by 


^  '  *n,p^^n!p*n.p!'\!p  in 


while  the  minimizing  forward  prediction  error  for  this 
selection  becomes 

-  *n.pt^n;p*n.pl''^n%lin  ^^^b) 

*^XY  in 

We  have  here  used  the  compact  matrix  product  repreeent- 
aClone 

Pw  ■  K  „  (13e) 

Xy  n.p'  n.p  n.p  n.p 


pC  ■  I  -  p 

*^XY  XY 

Since  we  are  only  Incereeced  In  the  optimal  ^  and 

,  wa  will  drop  cha  aymbol  and  aasua*  chat  x_ 

» n 

and  t*  ara  the  optimal  onaa  a*  given  by  aquation  (12) 

■  P  »  W 

Wa  may  aXao  daflna  cha  dalayad  backward  pradietion 
arror  vaecor  for  x^  by  ^ 

J  X  . 

■=p,n  -n  n,p  -p 

It  can  be  seen  chat  the  k^^  row  of  equation  (14) 
rcpreaents  a  prediction  of  x(k-p-l)  by  a  llnaar  combina¬ 
tion  of  cha  p  most  ioBodlaca  fuCura  valuaa  x(k-p) , 
x(k-p4-l) . x(k-l).Tha  raaulclng  error  In  this 

backward  prediction  la  d‘  (k) .  In  thla  caae  cha 
p,n 

optimum  vector  Is  Che  one  that  minimizes  cha 
quadratic  function 

8(i,)  ■  [d*  1+Wd*  (15) 

“1)  Pffl 

where  W  Is  defined  in  equation  (9) 

In  this  case  the  opcliwl  a  la  given  by 

a  *  -  -[Y  *  X  ^  (S'^^x..)  (16) 

P  n.p  n,p  n.p  — n 

In  a  similar  manner  to  Che  forward  prediction  error  case 
Che  optimal  estimate  of  sP^^Xn  is  specified  by 


-X  «- 

n.p  -P 


and  the  optimal  delayed  backward  prediction  error  vector 

X..5P+I  -g'^^..pC(sP*l  )  {17b) 

where  Pjy  and  given  by  equation  (13) .  Further¬ 

more,  It  can  be  shorn  that  ch*  optimal  ap  as  given  by 
equation  (16)  also  arises  by  approximating  p  Yule- 
Ualker  equations  In  a  asnner  similar  to  tha  approxlzM- 
clon  given  by  equation  (3)  for  ip. 

It  Is  clear  chat  the  forward  prediction  error  vector 

fx  and  the  aucoragreeeive  coefficient  vector  e.  are 
■^.n  — P 

Interchangeable  in  the  sense  chat  one  con  always  be 
found  from  the  ocher  using  equation  (7).  Similarly 
dx  and  Sp  are  Interchangeable  since  on*  can  always 
found  from  Che  ocher  using  equation  (14) .  It  la 
also  true  chat  the  2p  elemantt  of  ^  and  ^  sea  Inter¬ 
changeable  with  the  2p  elements  f*  ^(n) ,  f  J  |j(n) ..... 

f*  (n)  and  d*  (n) ,  d*  (n) ,  ....  d*  (n)  (1.*..  the 
pgO  IgO  ZgO  PgD 

nth  elements  of  chs  2p  prsdlccion  trror  vectors  ^ 

f*  . f*  and  d*  d?  . d*  ) . 

— 2,a  ’-^.n  — l.n  — 2,n  P,n 

Ue  will  show  this  last  fact  In  Section  VIII,  where 
we  will  also  see  chat  the  prediction  errors  lead  to  e 
lattice  filter  structure  which  la  related  to  the  auto¬ 
regressive  coefficient  vectors. 

In  Che  fast  cecucslve  algorithm,  the  sucoragtesalve 
coefficlanc  veccora  Op  and  ^  ara  not  directly  updated. 
Instead,  cha  prediction  error  elements  f?  (n) . 

1  gH 

f^  (n) ,  and  d?  (n) . d’^  (n)  are  updated.  Slnca 

p,n  p*o 

these  elements  are  Interchangeable  with  the  auto¬ 
regressive  coefficients,  there  la  no  Information  lost 
In  updating  only  the  prediction  error  elamaots.  How¬ 
ever,  cha  prediction  error  elements  may  be  updated  In 
a  computationally  efficient  manner,  requiring  0(p) 
mulcipllcaclona  and  additions  for  the  updata.  More¬ 
over,  as  wa  shall  see  later,  the  2p  prediction  arror 
elements  enable  us  to  f ind  all  of  the  ^  and  ^ 
vectors  for  ARMA  denominator  orders  from  I  to  p.  It 
is  for  these  reasons  that  we  choose  to  update  the 
‘^We  use  the  term  "delayed''  becaue*  although  tha  subacrlpt 
n  appears,  x(n)  is  never  used  in  (14).  Tha  undelayed  back¬ 
ward  prediction  error  vector  will  be  diacuseed  In  a  latar 
section. 
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pradlctlon  arror  elcaants. 

Vc  Bay  alao  obcaln  addlclonal  autoragraaalva  co- 
aff  iclaot  aatiaataa  slaiJLar  to  a  and  a  by  constdarlng 
tha  pradlctlon  arror  vaetora  asaoclacad  with  Cha  vaccor 
^  glvan  In  aquation  (4b).  Spaclflcally,  eba  opciaua 
forward  pradlctlon  arror  corraapondlng  to  2a  4a  daflaad 
aa 

•  V  4-  t  c*  (18a) 

-p,n  n,p  •=p 

uhara 

Sl  -  -tCn  (!«»»> 

ttfp  Q«P 

and  Xn,p  tn.p  and  /a  daflnad  In  aquation  (4).  Aa 
In  earlier  casea,  can  be  found  by  appronlaatlng  p 
Yule-Walker  eqtiatlona  or,  equivalently,  by  nlnlalalng 
tha  quadratic  functional 

Mcp)  -  (19) 

Correapondlng  to  thla  optlaal  autoragraaalva  co¬ 
efficient  vaccor  ^  wa  nay  alao  define  the  aaclatfta 

K'  -’^n.p  4  ■  "n,ptXn!p^a.pJ'\!p!(n 
-^X2n 

which.  In  cum,  glvea  rlae  to  cha  opclnua  error  vector 


where  and  P^^  are  defined  aa  in  equation  (13) 

Finally,  tha  opclnal  delayed  backward  prediction 
error  vector  for  2n  4*  daflnad  by 


d^  3*^^  +  Y  c 
2p,n  ■In  n,p 


V  •  -  lV,p’(n,pl‘\lp(s‘^ 


\gala,  w€  can  dafine  Che  predicted  value  of  X 

S%n*-  1fn.pC*t,Vn.P^'Vp<"‘^^i^ 

(22) 

and  It  followa  chat 

4..'- '« <“> 

Juac  aa  for  the  Xn  vector,  Che  2p  antltlas 

£{.„(") . ip,I(")’  =“ 

efficiently  updated  and  enable  ua  to  datemlna  the 
optloal  Cg  and  ^  coafflclenc  vaetora  for  all  AXKA 
nodal  denominator  orders  from  1  to  p.l 


III.  THE  HILBERT  SPACE  SETTING 

Tha  problem  of  recursively  updating  Che  prediction 
arror  vectors  in  Che  fast  algorithm  can  be  more  easily 
understood  by  casting  Che  problem  In  a  Hilbert  space 
setting.  Consider  the  o  dlswnslonal  coaiplax 
Euclidian  space 


H  •  (f' 


'  C  •  .  ■ .  *  C 


%rlth  Che  standard  vector  Inner  product  daflnad  by 
a 

+  r  * 

<x,2>  •  X  2  *  Z  *(4)  y(l) 
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’Throughout  the  remainder  of  the  paper,  the  symbol 
will  be  dropped  and  the  prediction  error  vectors 
fa,  da,  fP,  and  dP  are  assusied  to  be  cha  optimal  ones. 


We  note  that  cha  nsl  vaetora  yg.  and  ara 

all  aleasnts  of  S.  Moreover,  tha  p  columns  of  mstrls 
Xg  _  ara  also  elamencs  of  B.  The  sat  of  all  linear 
cottlnatlons  of  these  p  a laments  la  a  subapaea  of  B. 
idilch  we  denote  by  H^.  Similarly,  Hy  la  tha  subapaea 
spanned  by  the  p  columns  of  Tg^p. 

Let  us  now  consider  tha  forward  prediction  of  |^.  From 
aquation  (12a)  wa  saa  that  la  foensd  by  a  matrix 
multiplication  involving  ,Sa.  Tbs  matrix  Pn  la  saan  to 
be  a  linear  operator  on  tha  Hilbert  space  B.  It  la  ap¬ 
parent  that  Pjy  map*  alamants  of  B  Into  alamanta  In  tbs 
subapaea  Hg,  that  la 

'’XY=  «  *  “X  , 

Also,  It  is  evident  fron  aquation  (13)  that  Fgy  •  Pgy 
so  chat  the  operator  Fgy  la  a  projactlon  operator  onto 
tha  subspace  Mg.  In  general,  Pgg  Is  not  the  orthogonal 
projection  operator  onto  subapaea  Kg.  Inataad,  tha  as¬ 
sociated  direction  of  projection  Is  datarmlnad  by  tbs 
matrix  Yg^p.  It  can  ba  saan  from  aquation  (6)  that  tha 
dlrectlon’of  projection  of  Pgy  la  orthogonal  to  My. 
Thus,  P^  Is  cha  projection  operator  onto  the  subapaea 
Mg  along  My 4  (the  orthogonal  conplamant  of  My). 

With  chase  thoughts  In  adnd,  wa  can  provide  a  simple 
geonacrlc  Intarpracaclon  to  cha  four  error  vectors 
described  In  the  last  section.  In  particular,  the 
geonacrlc  relationship  beewaan  Xg,  and  ff  g  la 
depleted  In  Figure  1.  ’ 


Figure  1:  Geometric  Relationship  Betvean  Xg,  Xg,  and 
Che  optiswl  prediction  arror  ls~  ^ 

The  vector  ^  Is  seen  to  ba  that  projection  of  i^onco 
Mg  that  Is  orthogonal  to  My .  We  note  from  Figure  1  chat 

il,n  i  My 

or,  equivalently,  chat 


-P.n’  ® 


1.2 . . 


The  geometric  relationships  for  d^  g,  f2  g,  and  d£  gSre 
similar  to  Figure  1.  ’  ’  ’ 

Since  Pgy  ^  Pyg  ace  projection  operators,  so  are 
their  complamants  P^  and  P^.  It  follows  chat 

(P=  )^  ■  P=  C28a) 


We  can  also  see  from  equation  (13)  chat  ^XY  and  ^YX 
are  strongly  rslaced,  naiMly 

P„  -  [P„]^  (29a) 

and  P^  -  IP^jj)*  (29b) 

In  addition  to  the  four  pradlctlon  error  vectors, 
there  ara  four  inner  products  chat  ara  useful  In  deriving 
Che  fast  algorithm.  These  complex-valued  scalars  arc 
daflnad  as: 

“p.n  ^ 
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T  i  (d*  ] 

'p  n  ^  p ^ n * 


t«C 


IV.  THE  PROJECTION  OPERATOR  THEOREM 


1+0=  r«P+l 


(31) 

(32) 

(33) 


From  Che  results  of  the  last  section  It  la  apparanc 
chat  Che  various  prediction  error  vectors  and  scalars 
are  all  described  by  the  operators  Pj^y  and  Pyj^.  As  a 
new  data  point  x(n+l)  becomes  available,  we  desire  to 
update  Che  prediction  errors  vectors  and  scalars  In  a 
computationally  efficient  manner.  Because  the  opera¬ 
tors  Pj(y  and  ?Yx  used  repeatedly  and  their  struct¬ 
ures  change  as  new  data  points  become  available,  we 
prefer  to  update  them  and  Chen  obtain  updated  error 
vectors  by  applying  these  updated  projection  operators. 
Recursive  update  equations  for  Che  projection  operators 
Pyy  and  Pyx>  are  readily  obtained  by  appealing  Co  the 
following  theorem. 

Theorem  1.  (Projection  Operator  Theorem)  Let  A  and 
B  be  n«m  matrices.  Furthermore,  consider  Che  aug¬ 
mented  matrices  A  -  [Ala]  and  B  >  [B  I  In  which 
a  and  b  are  n^l  vectors.  If  [A  +Bl“^  and  [A  +B]“^ 
exist,  Chen  the  associated  projection  operator  corres¬ 


ponding  to  Che  augmented  matrices  Is  given  by 


(34) 

where 

Fab  -  A[B+A]'V 

(35a) 

-  ‘’ab 

(35b) 

_The  theorem  may  be  straightforwardly  proven  by  writing 
[A*B]  In  cerms  of  Schur  complements  and  performing 
some  matrix  algebra.  Alternatively,  the  theorem  may 
be  proven  using  Hilbert  space  concepts. 

Using  this  projection  operator  theorem  we  may  now 
obtain  all  the  necessary  equations  for  the  fast  recur¬ 
sive  algorithm.  Two  types  of  recursive  equations  are 
of  Interest.  First,  equations  are  needed  chat  provide 
the  m+l£t  order  prediction  error  vectors  In  terms  of 
Che  mill  order  errors.  These  are  called  order  update 
equations.  Second,  equations  are  needed  chat  enable 
us  to  update  the  prediction  errors  as  a  new  data  point 
becomes  available.  These  equations  are  referred  to  as 
the  time  update  recursions. 

These  two  secs  of  equations  are  derived  below. 


V.  ORDER  UPDATE  RECURSIONS 

In  this  section  the  order  update  equations  for  f^, 
fF,  1^,  qF,  u  and  v  are  derived  by  making  use  of  the 
projection  operator  theorem. 

Consider  first  the  error  vector  X’^oH-l  n  associated 
with  the  optimum  mS-lii.  order  autoregressive  co¬ 
efficients.  Here,  m  can  cake  on  any  value  In  the 
range  0,  I,  ....  p-1,  where  p  Is  the  desired  auto- 
rec'.esslve  coefficient  order.  From  equations  t^)  and 
Ulb).  we  see  that 


.X  _c 

t  ,  *  Pm  X 

(36a) 

where  X-X  ,  •  [X  1  s’^^’x  1 

n,ar»-l  n.io  .  -fl 

(36b) 

Y-Y  ,  ■  [Y  1  S®*^y  ] 

n.BM-l  n.m  .  •'■n 

(36c) 

Applvlng  the  projection  operator  theorem  to  (36a)  with 
A*X,  B-Y,  a  -  S®+lx,  and  b-S^t+l-y,  we  have 

^rwl.n  ^XY  - 


-  P 


(S°^’'x)((S°^‘y)"p^(s'^^x)l  5. 


’m,n 


m,n 


‘XY 
(37) 


As  mentioned  earlier,  to  Implement  the  recursive  al¬ 
gorithm  we  only  need  the  error  element  at  tisM  n,  that 

Is  f^,  (n) .  From  equation  (37)  we  see  that 

mri , n  . 


irt-l,n 


(n) 


f'^  (n)  -  d*  (n) 

m  n  V"  m  n 

m,n 


(38) 


Equation  (38)  la  the  desired  order  update  equation 
for  fx. 

In  a  similar  manner,  the  order  update  equation  for 
fP  Is  found  to  be 


•^iirt-l,n 


’’7x2n 


where  X  and  Y  are  defined  In  equation  (36). 
the  projection  operator  theorem  yields 


-m+l ,  n  -m , 


+m.n 

V 

m,n 


dP 


(39) 
Applying 

(40) 


The  niil  component  of  equation  (40)  Is  the  desired 
order  update  equation  for  the  forward  y  prediction 
error,  that  is 

01  pH  ''01  pH  IllpD 


The  order  update  equations  for  the  delayed  backward 
prediction  error  vectors  may  be  similarly  derived. 

These  equations  are,  however,  not  as  useful  as  the 
combined  order  and  time  update  equations.  The  combined 
update  equations  give  d*  ,  and  dP  ,  In  terms  of 

tBri  j  nri  QrX  p  Hri 

X  y 

d  and  d  In  deriving  the  combined  order  and  time 

-HOpn  ^Bpii  * 

update  equation  for  the  delayed  backward  ^  estimate, 
we  first  note  that  from  aquation  (17b) 


jX 

^m+l,n+l 

(42a) 

where 

r.a » 

X  ■ 

n+1, 

s 

.m+I 

1  X 

i 

(42b) 

L,  Q>®* 

-“j 

1 

Y  -  P-., 

m  I 

(42c) 

'n,ml 


We  now  apply  the  proj^ection  operator  theorem  to 
equation  (42a)  with  A>X  and  B>Y  In  equation  (34). 
After  some  simple  algebraic  manipulation,  we  get 


-m+1 ,  n+1 


r..o..' 

[^.n 


(43) 


m,n 


The  irflat  component  of  equation  (43)  yields  the  desired 
update  equation 


Vl,rH.l^'^l^ 


Sp  Q 


n,n  -X 
*  ^ 


n»n 


(n) 


(44) 


m,n 


The  delayed  backward  time  and  order  update  equation 
for  ^  Is  derived  In  a  similar  manner.  The  details  are 


omitted,  but  It  Is  readily  shown  that 

Cl.n+l^'**^^  ■  '‘m,n<"> 

-  -21^  fP  (n) 
u  ^  nipn 

(45) 

(Dpn  * 

Finally,  the  order  update 

equations  for  the 

scalars 

Um,n  s“d  v„  „  are  derived. 

Fron  equation  (32) 

“n+Un  ■^n'’5S 

(46a) 

wh«r«  X  ■  X  .  ■  (X 

n,nH*l  n,m  . 

(S“^^x  )] 

(46b) 

y  -  Y  -  [Y 

0,111+1  n,m  . 

(46c) 
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Applying  the  projection  operator  theoraa  glvaa  Che 
daalred  result 

“irt.l.n  •  i 

-i+pC  (S**'^*)  (S**’^y)^p'  X 


u  .  •  u 

*fl.n  a.n 


Also,  from  equation  (33) 


Vl.cH-l  ■  <^®> 

where  X  and  ¥  are  defined  In  equation  (46) .  Applying 
the  projection  operator  theorea  yields  Che  coablned 
order  and  tlaa  update  equation  for  v. 


m+l,n4-l  a,n 


VI.  TIKE  UPDATE  EQUATIOHS 

The  reaalnlng  recursive  equations  update  the  forward 
and  delayed  backward  errors  as  a  new  data  point  Is 
obtained.  For  this  reason  these  equations  are  called 
clae  update  equations. 

When  a  new  data  point  becoaes  available,  Che  effect 
on  Che  prediction  error  vectors  is  to  append  a  row  to 
Che  boccoa  of  their  defining  aacrlx  equation  [aae,  for 
exaaple,  equation  (7)].  Appending  a  row  to  Che 
aacrlces  Xq  n  and  Y„^o,  does  not  seea  to  fit  In  Che 
fraaework  o^  Che  proj’cctlon  operator  cheorea.  In  which 
coluana  are  appended  to  Xq^o  And  Ya,a|.  It  turns  out, 
however,  that  we  can  accoapllah  the  cask  of  annihilat¬ 
ing  a  row  In  Che  error  vector  aacrlx  equation  by 
appending  to  Xq^,,  and  Y^^,  the  qCh  basis  vector  defined 

by  ’  I 

e^  -  (0  0  ...  0  11  (30) 

To  see  how  this  works,  let  us  consider  as  an  exaaple 
the  forward  prediction  error  vector  for  X;n.  If  we 
append  the  nxl  vector  ^  Co  Xn^g,  In  aquation  (7)  we 
have 


f*  (n)l 
Lm,n  j 


x(n-l)  x(n-a) 


a.n  •  Sn 


Equation  (33b)  Is  satisfied  only  if  we  force 
(n)  •  0.  This  can  always  be  done  because  Che 

o«  o 

scalar  C  appears  only  in  the  last  row  of  equation  (51). 


la  parClculaFp  -  0  If  we  choose 

a 

e  -  -  Z  a_(l)x(n-l)  (54) 

1-1 

Since  {(.nCh)  -0,  equation  (53a)  is  seen  Co  depend 
only  on  the  first  n-1  coagKiaencs  of  the  vectors.  It  is 
easily  seen,  then,  chat  (k)  for  k-1,2 . n-1  are 

deterained  using  the  first  n-1  rows  nf  (51)  in  such  a 
aanner  that  the  vector  (f*  _(!),. ..,ff  _(n-l)j'  is 

orthogonal  to  the  first  n-1  coaponencs  of  each  colusB 
of  Yn,B.  But  this  Is  exactly  the  problesi  of  decenin- 
Ing  Che  forward  prediction  error  based  on  n-l  data 
points.  Thus,  we  see  chat 


similar  argioMncs  show  that  this  time  annihilation 
property  alao  holds  for  n>  IJ  q  with  re¬ 

sulting  formulas  similar  to  equation  (35).  We  finally 
note  that  the  scalars  o^.n.  U*  n«  sre 

formed  as  Inner  products  of  the  prediction  error 
vectors.  Since  the  last  element  of  t*  (or  f^  d*  , 

d^  )  is  zero,  it  follows  that 


Vn-l  ■ 


and  similarly  for  T„  „,  ,  and 

With  these  thoughts  In  mind  we  are  In  a  position  to 
derive  time  update  equations  for  a,  t,  u.  and  v.  First, 
let  us  define  the  augsMnted  matrices 

*-I*n.«:5nl  »7a) 

(57b) 

Then  It  follows  chat 

°m.n-l  ■  (^8) 

where  X  and  f  are  defined  by  equation  (57).  Application 
of  Che  projection  operator  theorem  yields 

°m,n-l  “  ^[’'xY  ■  ^XY 

'’m.n-l  -  “m.n  ’  ^m.n^OX^  -^m,nJ''-*m.n(“>  (^9) 


where  1  -  y_  _  4  •- 

m,n  -  — n  XY  —n 


By  rewriting  equation  (59)  we  arrive  at  Che  desired 
time  update  equation 

[f*  (n)!^  „(n)l  (61) 

BgP _ BaO _ 


where  ^  is  used  instead  of  ^  ro  indicate  Che  presence 
of  Che  e„  vector.  The  optimal  f^^n  vector  la  given  by 

^  n  ‘  ‘’K  in 
where  ^ 

^  ;  Sn' 

<520 

From  equation  (27b)  we  know  chat 

'  fm.n  •  >  -  0  •  i-1.2 . »  (53a) 

<  f2  n  .  e  >  -  0  (53b) 


0  -0  ,  +• 
m,n  m,n-l 


The  time  update  equation  for  t  Is  found  In  a  similar 
manner  by  using  2  for  A,  ?  for  B,  S^A^^j^for  a,  and  ^ 
for  b  In  equation  (34)  to  yield 

n(“>l  (“> 

.  .  ,  +  _n,n _ mjn _ 

m,n  m,n-l  1  "  Y_  - 

B«n 

ThB  updatB  •quBtlons  for  u  and  v  art  found  to  bt 


U  •  U  1 

m.n  B,n<-I 


V  —  V  , 

BpO  Bpn-1 


(f*  „(n)l*[f.  „(“)) 

BsO  _ BtP 


(d*  „(n)l*[d^  „(n)] 

Bp  P _ BbP 
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by  using  arguments  anelogous  to  those  used  In  deriving 
equations  (61)  and  (62). 

Finally,  update  equations  for  Ygi,n  needed.  To 
obtain  an  order  update  aquation,  we  note  that 

1  *  Y  *  e  e 

’'m+l,n  =n  XY  ®n 


(65a) 


where 


(65b) 

(65c) 


Applying  the  projection  operator  theorem, 
1-Y  "e*p'^  e 

’nri-l,n  ^ 


1  -  Y, 


m,n 


[<i:  An)r[d^  ^(n)] 


m,n 


or«  equlvalencly 


(d^  „(n)l  [d^„(n)] 
De  n  n. n 

Y  a  Y  ^  ■!  fa  I  III  1 

oH-l ,  n  m ,  n  * 

u 

m,  n 


(66) 


(67) 


The  (undelayed)  backward  errors  are  defined  by 


b"’ 

-P,n 


-p,n 


S^x  +  [S°x  :  S^x 


P-1 


:  s  ; 


:  s’ 

;  sP-^ 


(71) 

(72) 


where  and  c^  are  p><i  aucoregrcaslve  coefficient 
vectors;  By  comparing  equations  (71)  and  (72)  with 
equations  (14)  and  (21),  It  la  readily  seen  that  the 
(optimal)  backward  and  delayed  backward  error  vectors  are 
related  by  r- 


and 


d^^  - 


b’‘  , 

— tn,n-i 


.p... 


o,n-l 


(73a) 


(73b) 


The  nfh  components  of  equations  (73)  yield  the 
desired  results 


:D«n 


,  (n-l) 


m,n-l 

d^  (n)  -  b^  , (n-l) 
01  n  i&  ri*l  ' 


(74a) 

(74b) 


To  obtain  the  tine  and  order  update  equation  for  Y 
we  note  chat 


^  -  '^m* 
where 

X  •  .X 


•  e 


n+1 ,  oM-l 


-n+1 

(68a) 

] 

■*  1 

(68b) 

n,m_| 

'  '1 

1 

(68c) 

n,mj 

_ro,  j  i' 

^n+l,iB+l  ’  '.  Y 


Since  Che  first  element  of  e  .  Is  aero.  It  follows 
chat 

^  -  "m^l,n+l  '  •£n‘'’lY 


where 


X  -  [X  ;  X  ] 
n,m.— n 


(69b) 

Applying  the  projection  operator  theorem  yields 


1-Y, 


nrt-1 ,  n+1 


;  _  Y  .  _ 1°,” 

'm,n  * 


^m+l,n+l  ^m,n  ^ 


n,n 


K  n^nJnf'  „(n)l 

01  s  n  B  t  ^ 


(70) 


Thus,  using  equation  (74)  we  can  replace  the  delayed 
backward  error  terms  In  Che  recursive  equations  by 
backward  error  terms. 

It  Is  also  helpful  to  replace  v  by  another  scalar 
defined  by 


0  -  [b*  )*(b^  ) 

m,n  -m,a  -m,n 


(75) 


By  considering  Che  defining  equation  for  v,  equation 
(33),  along  with  equation  (73)  It  Is  a  simple  matter 
CO  show  that 


ID  ,  •  V 

a, n-l  a^n 


(76) 


Table  1  summarizes  Che  complete  set  of  update  equations 
with  d*,  d7,  and  u  replaced  by  b*,  bY,  and  u. 

The  Initial  conditions  for  the  various  quantities 
In  Table  1  are  obtained  by  considering  the  defining 
equations  for  these  quantities.  For  exasqtla,  from 

equation  (24)  we  see  chat  If  m  ■  0  then  f!  ■  x  or 

■^*11 


fp_^(n)  -  x(n) 

From  equation  (71),  It  also  follows  that 

b*_n(n)  ■  x(n) 

Similarly,  we  can  show  that 
fg  jj(n)  ■  y(n)  ■  x(n-q) 

bn  _(n)  •  y(n)  •  x(n-q) 

U|  n 


(77) 

(78) 

(79) 

(80) 


These  six  time  update  equations  (61),  (62),  (63),  (64), 
(67),  and  (70)  along  with  the  six  order  update 
equations  (38),  (41),  (44),  (45),  (47)  and  (49)  obtain¬ 
ed  In  Che  last  section  comprise  the  fast  recursive 
algorithm. 


VII.  SUMMARY  OF  THE  ALGORITHM 

AC  this  point  all  necessary  equations  for  Implement¬ 
ing  Che  fast  recursive  algorithm  have  been  derived.  In 
this  section  Initial  conditions  are  discussed  and  a 
procedure  for  Implementing  the  recursive  equations  Is 
given. 

The  implementation  of  Che  algorithm  Is  conceptually 
easier  If  we  replace  the  delayed  backward  error 
quantities  by  (undelayed)  backward  error  quantities. 


Furthermore,  since  vn  is  the  zero  vector  for  n^q. 
It  follows  from  equations  (30)  -  (33)  that 


0,n 


T„  •  Un  *  wn  “0  for  n  <  q 

Uftl  U|Q  UpD  ~ 


(81) 


It  Is  clear  from  Table  1  chat  no  order  update 
equations  exist  from  o  and  t.  Therefore,  Initial 
conditions  are  needed  for  Op,n  Tp,Q  for  each  p  and 
some  corresponding  n.  A  little  thought  will  convince 
one  that  f*  q^pfq+f)  •  0  by  using  an  argument  similar 

to  Che  time  annihilation  argument  presented  In  Section 
VI.  Moreover,  we  can  also  show  chat  d’^  ^  (q+p)  ■  0. 

P,q+1> 

by  using  this  type  of  argument.  Thus, 


for  m  •  0,1, . . . ,p 


(82) 
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(T-1) 


fm,n<“> 

- 

b*  ,(n-l)  a*  /u*  , 

(T-1) 

- 

(T-2) 

‘>^.n-l(“-^> 

- 

f“  (n)  t*  /U* 

Syfi  fQyQ  Byll 

(T-3) 

‘>1  n 

iDg  n**i 

- 

f^  (n)  0  /u 

BfH  Byd  Byd 

(T-4) 

u 

fflgii 

- 

0  T  /u  , 

Byd  Byd  Byd^l 

(T-5) 

- 

3  T  /u 

Byd  Byd  Byd 

(T-6) 

a 

B,n 

+• 

Byd  Byd**l  Byd 

(T-7) 

T 

m,n 

^m,n-l 

+ 

Byd^X  Byd  Byd 

(T-8) 

**31,  n 

'^m,n-l 

tf*  „(n)]*[f^  ^(n)I/(l->  ) 

Byd  By  d  Byd 

(T-9) 

j 

m,n 

“'m,n-l 

+ 

tb*  (n)]*[b^  „(n))/(l  -Y„  „) 

Byd  Byd  Byd 

(T-10) 

^awl,n 

'^m,n 

+ 

Byd^l  Byd*i  Byd"! 

(T-11) 

'^orFl ,  n+1 

'^m,n 

(f*  „ 

Byd  By  d  B,  d 

(T-12) 

where  p  is  Che  desired  (oaxlmun)  autoregressive  co¬ 
efficient  order. 

Finally,  when  m  •  0, 

^O.n  -  °  («) 

Although  ocher  initial  conditions  nay  be  obtained, 
these  initial  conditions  are  Che  only  ones  needed  to 
ispleoent  the  algorlchn. 

The  Implementation  of  Che  update  formulas  can  be 
divided  into  three  parts.  First,  for  n^q  the  vector 
^  is  the  zero  vector,  so  no  operations  are  performed. 
For  q'*.l^n;^q.*-p.»'l,  Che  maximum  order  m  Chat  can  be 
used  Is  n-q  -1.  In  this  time  interval,  as  a  new  data 
point  arrives  not  only  are  time  updates  performed,  but 
also  Che  model  order  is  increased.  For  n>q')-p-*-l,  the 
model  order  remains  at  p  and  the  time  updates  only  are 
performed . 

The  implementation  of  the  algorithm  for  n^q-t-1  is 
summarized  below. 

As  the  new  data  point  becomes  avaiiable; 

1)  Set  n  -  n+1 

2)  These  quantities  are  available  from  the  last 
Iteration: 

f*  ,i:n-l),f^  „  .(n-D'bl  „  ,("-!) 

tUsn—i  Qifi'*!  ni,n»l 


Table  1:  Summary  of  Update  Equations 

4)  For  each  m  ■  0,1, . . . ,min[p-l,n-q-31  find 


'm,n-l'  'Ti,n-1’ “m,n-l'‘^m,n-l 


3)  New  initial  conditions: 


If  n^p+q  +  2. 


'0,n‘"> 

‘>S,n^"> 

^O.n  ‘ 

'n-q-l,n-l 


x(n) 

x(n) 

x(n-q) 

x(n-q) 


T  ,  ,  •  0 

n-q-l,n-l 


o  ,T  ,w  ,u>  using  (T-7)-(T-10) 

OpQ  nipQ  EDiU  o^n  ” 

(T-l)-(T-4) 

Vl,n  4 

5)  For  m  •  min[p,n-q-2]  find 

o  ,T  ,u  ,u  using  (T-7)-(T-10) 

tBpIV  IBpD  OpD  Bpll  * 

6)  If  n  <  p  s-q  + 1  we  need  to  add  a  filter  order. 
Set  ra  -  n-q-2. 

Find: 

Cl,n(")>Cl.n(">-Cl,n(">>'’L.l,n("> 

(T-l)-{T-4) 

W.n 

“nrt-l.n’-'oH.l.n  •  <^-6) 


for 

Set  m 

a“0,l 

Find: 

•  •  •  y 

0 

min(p. 

m,n 

n-q-2) 

At  ch: 

0  ,  ■ T  ,"0 

m,n-l  m,n-l 


to  arrive. 

It  is  clear  from  Che  above  suanary  that  0(p)  multi¬ 
plications  and  additions  are  required  to  update  the 
prediction  errors.  More  specifically,  in  Che  time  up¬ 
date  mode  [i.e.,  when  nsps-q  +  l  so  no  filter  orders 
need  to  be  added]  14p  multiplications  and  lOp  additions 
are  performed  per  update.  In  the  tisn  and  order  update 
mode,  [i.e.,  when  n^ps-q+l),  17p  multiplications  and 
13p  additions  are  performed.  It  should  be  noted  chat 
this  computational  requirement  may  be  significantly 
reduced  by  using  a  normalized  lattice  form  similar  to 
chat  in  [2]. 
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\ail.  THE  LATTICE  STRUCTUBE 


The  fast  recucalve  algorlcha  lends  itself  to  a  digi¬ 
tal  filter  structure  known  ss  the  lattice  filter.  This 
structure  can  he  aeon  by  considering  equstlons  (T-1)  - 
(T-A).  It  Is  clearly  seen  that  the  filter  that  pro¬ 
duces  the  four  outputa  f*  _(n),  (n)  and 


a.n  a,n 

(n)  is  given  In  Figure  2,  where 
a  n.o 


y  <a«o 

a-'  •  — *— 


Figure  2: 

Filter  Realization  of  Error  Order  Update  Equations 

Thus,  the  entire  plh  order  filter  Is  given  by  p  stages 
of  filter  shown  In  Figure  2.  This  Is  depicted  In 
Figure  3 


p  stages 


Figure  3: 

Realization  of  the  p^*'  Order  Filter 

In  Che  startup  phase  of  the  algorithn  [l.e.,  when 
the  filter  order  Is  less  chan  Che  desired  order]  Che 


filter  of  Figure  3  begins  with  one  stage,  and  success¬ 
ive  stages  are  added  as  new  data  points  arrive.  Alter¬ 
nately,  all  p  stages  nay  be  la  place  at  the  beginning: 
however,  the  four  nulclpller  sections  of  each  unused 
stage  is  set  to  zero  until  chat  stage  Is  to  be  used. 

The  lattice  filter  structure  at  Figure  3  nicely 
depicts  Che  relationship  between  the  autoregressive 
coefficients  and  Che  prediction  errors.  To  see  this 
relationship,  let  us  denote  the  transfer  functions  fron 
x(n)  to  fi^  (o)  and  from  z(n)  to  b>  (n)  by  F*(x)  and 

■pH  ■pH  ■ 

Bj|(z),  respectively.  Let  us  also  define  Che  auto¬ 
regressive  coefficient  transfer  functions 

A  (r)  -  1  +  a  (l)z“^  +. ..+s  (■)z"*  (88) 

B  ■  ■ 

i  (z)  -  a  (0)  +  a  (l)z"^  +...+S  (ni-l)z"^^+z"" 

■  B  B  B 

where  Che  84  and  a^^  coefficients  are  defined  by 
equations  (5)  and  (71)  respectively.  A  little  thought 
will  convince  one  chat  the  transfer  functions  are 
related  by 

F*(z)  -  A.(z)  (90) 

B  B 


B*(z)  -  A  (z) 
a  B 


(91) 


With  this  thought  In  nlnd  we  are  now  able  to  derive 
equations  Chat  relate  the  error  elements  to  Che  auto¬ 
regressive  coefficients.  From  Figure  2  it  Is  clear 
chat 

C  -  Cl  n(“)  n**™  1  „  1  <”) 

BpD  a^lpH  B^ipXi  a**Xpii**l 

BpO  Oh^ipa  IIt**lpD  B^lpfl^I 

By  caking  z  cransforos  of  equations  (92)  and  (93)  and 
using  equations  (90)  and  (91)  we  get 

•  Vl<*>  ^*'^Cl.nVl(*> 

A^(z)  -  *'^Cl.nVl^*>  Vl^*^ 

Furthermore,  we  can  see  from  Figure  3  chat 

Ap(z)  -  Ap(z)  -  1  (96) 

and  from  equations  (94)  and  (95)  chat 

•  ®n,n  (97) 

n 

a  Bp  Q 

Given  the  error  elements  1*^  (n)  and  b*  (n)  for 

apEi  Bpii 

m»I,2,...,p  we  can  use  equations  (92)  and  (93)  to  get 
Che  o2,n  eod  n  elements,  then  use  equations  (94)-(9S) 
CO  get  Che  a,(k)  and  ag,(k)  elements  for  m>l,2,...,p 
and  k  •  0,1,2, ... ,n.  Thus,  we  are  able  to  obtain  Che 
autoregressive  coefficients  for  all  model  orders  from 
1  to  p. 

Similarly,  given  Che  ap(k)  and  a-fk)  coefficients  for 
k-0,l,...,p,  we  can  use  (94),  (95),  (97),  and  (98)  to 
get  the  Op,n  ®P,n  “*•*  ^**e”  P  "1*^  order  auco- 
regrasslve’coeff iciants  by  working  down  from  order  p 
to  order  0.  Then,  using  (92)  and  (93)  %ra  can  obtain 
the  desired  error  elements. 

Thus,  using  equations  (92)-(98)  we  arc  able  to  con¬ 
vert  back  and  forth  between  the  lattice  error  elements 
and  the  optimal  autoregressive  coefficients. 

As  a  final  note,  when  q*0  (l.e.,  when  an  auto¬ 
regressive  model  Is  chosen)  It  can  be  seen  In  Figure  3 
chat  Che  cop  half  and  bottom  half  of  the  filter  are 
equivalent.  In  this  case,  only  one  half  la  needed,  and 
the  filter  In  Figure  3  degenerates  to  the  AR  lattice 
filter  described  In  [1]. 
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IX.  CONCLUSIONS 

In  chls  paper  we  have  presented  a  recursive  algorithm 
for  obtaining  the  autoregressive  coefficients  of  an 
AJtMA  model.  The  recursive  algorithm  la  based  on  the 
prewindowed  version  of  the  high  performance  method  of 
ARNA  spectral  estimation  as  described  In  Part  1.  The 
recursive  algorithm  Is  computationally  fast,  requiring 
0(p)  additions  and  multiplications  to  update  the  para¬ 
meters.  Moreover,  the  algorithm  can  be  Implsmented 
using  a  lattice  filter  structure  offering  numerical 
robustness  and  nice  convergence  properties  associated 
with  lattice  type  algorithms. 

We  have  not  yet  discussed  the  problem  of  recursively 
estimating  the  moving  average  coefficients  In  the  ABMA 
modal.  Ua  do  not  at  this  time  have  such  an  algorithm. 
However,  it  Is  worth  noting  chat  Che  moving  average 
Information  Is  present  In  the  output  prediction  error 
sequences,  and  the  utilization  of  this  Information  for 
moving  average  coefficient  estimation  la  currently 
under  study.  Another  area  under  study  la  the  use  of 
various  normalization  procedures  to  effect  a  decrease 
In  computational  requirements  and  In  sensitivity. 

Finally,  we  note  chat  the  recursive  algorithm  pre¬ 
sented  here  Is  based  on  approximating  a  sec  of  p  Yule- 
Walker  equations.  It  has  been  shown  In  Part  1  of  this 
paper  that  for  short  data  lengths.  Improved  spectral 
estimates  result  from  using  more  than  p  Yule-Walker 
equations.  For  chose  cases  In  which  the  amount  of  data 
Is  small,  a  fast  recursive  algorithm  baaed  on  the 
jpproximacion  of  c  >p  Yule-Walker  equations  would  often 
prove  useful.  Such  an  algorithm  Is  currently  being 
pursued . 
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